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ABSTRACT. We compare the numerical solutions of the 2+1 equivariant Wave Map problem computed with 
the symplectic, constraint respecting Rattle algorithm and the well known fourth order Runge-Kutta method. 
We show the advantages of the Rattle algorithm for constrained system compared to the free evolution with 
the Runge-Kutta method. We also present an expression, which represents the energy loss due to constraint 
violation. Taking this expression into account we can achieve energy conservation for the Runge-Kutta scheme, 
which is better than with the Rattle method. Using the symplectic scheme with constraint enforcement we can 
reproduce previous calculations of the equivariant case without imposing the symmetry explicitly, thereby 
confirming that the critical behaviour is stable. 



1. Introduction 

Down to the present day, the use of symplectic, or more general geometric, numerical integration methods 
in the field of (general) relativistic (field) equations is not very common. To our knowledge the first work 
was done in [4] and [3] to solve ODEs in the context of cosmological space-times. The step to field 
equations was done in [27] and [28]. A symplectic integrator for post Newtonian equations was developed 
only recently [ ]. A reason why this has not been done in the past very often, is the different symplectic 
structure of General Relativity compared to usual Hamiltonian systems [ ]. 

Here, we want to present another contribution to the problem of solving relativistic field equations with 
the help of symplectic integrators. We compare the numerical solution of a 2+1 dimensional Wave Map 
system solved with a symplectic, constraint preserving (at least in a given accuracy) integration scheme 
with the standard fourth order Runge-Kutta method. We will show that we can reproduce the qualitative 
results gained by Bizoh et.al. [ ] and Isenberg and Liebling [ ] in the so called equivariant case. 

The outline of this article is the following: in section 2 we give a brief description of the symplectic 
integration method, which we will use for our numerical studies. As already mentioned, we solve a Wave 
Map system to show the differences between the integration methods. So in section 3 we describe Wave 
Maps in general and concentrate then on the 2+1 equivariant case. In section 4 we derive a correction 
term for the energy, which is caused by the constraint violation during the numerical time evolution. The 
numerical setup, including the choice of the initial data and the boundary conditions are presented in section 
5. Finally, in section 6, we present some results of our numerical studies. 

Our numerical results can be put into three independent categories: 

(a) The comparison of the standard fourth order Runge-Kutta method with the Rattle method. 

(b) The correction of the energy by taking the constraint violation into account. 

(c) The blow up dynamics of the 2+1 equivariant Wave Map, where we reproduce known results. 

2. The Rattle Method 

The Rattle method is a numerical integration method for constrained ODE systems. It is a further develop- 
ment of the Shake algorithm [ ] (which itself is based on the well known Stormer-Verlet method [ ]) to 
take the presence of constraints into account. It was developed by Andersen [ ] in the context of molecular 
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dynamics. Later, it was shown by Leimkuhler and Skeel [ ] that the Rattle method belongs to the group 
of the symplectic integrators. 

Symplectic integrators are a certain class of geometric integrators [9, 1 0, 2 1 , 24] . Their characteristic feature 
is the conservation of a geometric structure or quantity. Such quantities can be for example the energy, the 
momentum, the phase space volume or, like in the case of the symplectic integrators, the symplectic 2-form 
co = dq a Adp a . 

In many branches of physics and applied mathematics, including celestial mechanics, molecular dynamics 
and quantum mechanics symplectic integrators have been used very successfully for years. Very often 
the numerical results, especially in terms of long-term energy and (angular) momentum conservation, are 
much better compared to non-symplectic integrators [10]. 

Now we give a brief description of the Rattle method. For further details we refer to the original paper [ ]. 
We consider a Lagrangian system with holonomic constraints § A and an action of the following form 1 : 

*t[q,4,X] = J (^M kj q k qj-V(q)+X A <l> A (q))dt 

which leads to the Euler-Lagrange equations 

j t (M Jk q k )+djV-X A dj<j> A =0, 

(1) <p A (q)=0. 

The constraint equation (1) implies also an additional constraint on the velocities: 

W A (q,q):=^<j> A (q)=qJd j (l> A = 0. 

With M we denote the mass matrix, with q the generalized coordinates, q their corresponding velocities, Xa 
are the Lagrangian multipliers and V is the potential. With the small Latin indices we label the dynamical 
variables and with the capital indices the constraints. We use here and throughout the whole article the 
common dot notation for time derivatives: / := d t f. 

A time-step with the Rattle method starts with a Stormer-Verlet step using the unconstrained equations. In 
general the result does not satisfy the constraints. Therefore, the next step is to compute the appropriate 
amount of the constraint force — Xa djQ A , which is necessary to push the estimated solution to the constraint 
surface. This is done by first correcting the coordinates qi by iteratively solving the constraint <j> = to 
determine the values of the Lagrange multipliers. Once these constraints are satisfied to a given accuracy 
the velocity constraints y A = are solved in a similar way to correct the velocities. 

3. Wave Maps 

Wave Maps are systems of generalized wave equations. They appear in several physical models, for exam- 
ple as nonlinear sigma models in particle physics [12] as well as the Einstein equations for a certain class 
of electro- vacuum space-times [2] or in cosmological models [ ]. On the other hand, Wave Maps are 
of mathematical interest, because of their nonlinear structure and the possibility for singularity formation. 
Wave Maps are also the simplest example for geometric wave equations [32]. 

3.1. General. In general, a Wave Map describes a mapping from the m + 1 -dimensional Lorentzian man- 
ifold j% into a n-dimensional Riemannian manifold jY . Usually we distinguish between two different 
kinds of formulations for the Wave Maps. The first, so-called intrinsic formulation, describes the mapping 
directly from the base manifold into the target manifold JV . 

© : Jl ->• JV 

X^r^= &(x) . 



It is only possible for systems with holonomic constraints to write the action in this form. 
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We denote here and in the following the coordinates in j$ withx = (x°,x , .. and with | = (| ...,§") 
the coordinates in jV . Further, g ab is the metric on j& and hpy the metric on jV , so small Latin indices 
refer to manifold ^# and small Greek indices to jV . 

The action for the intrinsic formulation of the Wave Map is 

^[0,V0]=/ g ab V a &Pv b &rhp r ^/\g\d m + l x. 
We denote with V a the covariant derivative with respect to the metric g a b on and with g the determinant 

of gab- 

The variation of the action with respect to the field variables 0^ gives us the equations of motion: 

(2) n g &P + r^ MV v a @v V„0 V = . 

With □„ = V a V a we denote the d'Alembert operator (wave operator) on the base manifold and with 
r^ v = r^„ v (0(x)) the Christoffel symbols corresponding to the metric h^y on jV. Equation (2) shows 
clearly that Wave Maps are generalizations of the wave equation (jV = R) on the one side and the geodesic 
equation (ytt = K.) on the other side. 

A different way to describe Wave Maps is the so-called extrinsic formulation. In this formulation we 
assume the target manifold jV isometrically embedded into a higher dimensional euclidean space W , 
p > n. 

U: Jt^jV^W 
X H> z = U{x) 

with z = z l , ■■■,z p the coordinates on R p and the notation for the variables in ^# being the same as for the 
intrinsic formulation. In this case we need additional conditions <p (z) = 0, which ensure that the image of 
the mapping lies in the sub-manifold jV of W' '. We can do that by attaching the constraint expressions 
(j)(z) to the action via Lagrangian multipliers. For the following purpose we can assume that it is enough 
that p = n+ 1, which means that we need only one equation to describe the embedding of jY into R" + . 
We write for the action of the extrinsic formulation 

(3) */[U, VU,X] = J (g ab V a U A V b U B 8 AB + 2A0) y/\f\ d" t+1 x . 

The capital Latin indices refer to and 8ab denotes the Euclidean metric on K n+ . 
Again, extremization of this action results in the equations of motion for this formulation: 

a g u A +Xd A (j) = o, 
<Kz) = o. 

3.2. 2+1 Wave Map. In the following we will concentrate on a Wave Map with the base manifold ^ — 
M 2+1 , the 2+1 dimensional Minkowski space-time and the target manifold J/ = § 2 , the 2-sphere. 

Intrinsic Formulation. We will use the intrinsic formulation to discuss the basic properties of the 2+1 
Wave Map. In this formulation, we write for the map: 

: M 2+1 =R 2 xR^§ 2 

{x°,x\x 2 )^ {■&(x°,x\x 2 ),(p{x Q ,x\x 2 )) . 

With the common notation for spherical coordinates on the 2-sphere we write for the line element 

ds 2 2 =d# 2 +sin 2 #d<p 2 . 

The equations of motion take the following form: 

(4) □ t! #-sin#cos#V>V fl (p = 0, 

(5) □ ? (p + 2coti?V < 'jJV„9 =0. 

Further, we choose polar coordinates on the Minkowski space-time with the metric 

di 2 2+1 = g ab dx a dx h = At 2 - dr 2 - r 2 da 2 
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where r is the radial and a the angular coordinate. The most commonly studied situation for this Wave 
Map is the so-called equivariant case. In this case one sets #(f , r, <r) = #(f , r) and <p(f , r, a) = a. The latter 
assumption means that the Wave Map maps rotations on Minkowski space around the origin to rotations of 
the sphere around the 3-axis with the same angle. In this case, the equation (5) is identically satisfied and 
(4) reduces to 

(6) #_^_V + *W =0 . 

r Zr A 

By confining to the equivariant case one needs to solve only one equation, but at the cost of having to deal 
with a coordinate singularity at r = 0. This equation together with the initial data 

(7) #(0,r) = tfo(r) and #(0,r) = ^(r) 

is the Cauchy problem, which was studied in [ ]. Equation (6) has two non-trivial, static solutions with 

0(r) G [0,71]: 

(8) # s (r) = 2arctan(r ±1 ). 

The local well posedness of the Cauchy problem for Wave Maps with general target manifold, was proven 
in [13] and [15]. The study of the global well posedness started with the investigation of the equivariant 
case. This was done in [6, 7, 34-36]. The problem of global well-posedness for the non-symmetric 2+1 
Wave Map with small initial data 2 was first addressed by Tataru [42]. Later work include [14,33,40,41,44]. 
For the case of a hyperbolic target manifold, see [16, 17]. Just recently the global well-posedness for the 
Cauchy problem in 2+1 dimensions for large energy initial data was first proven by Krieger and Schlag 
for the hyperbolic target [ ] and by Sterbenz and Tataru for the § 2 case [38, 39]. Finally, we want to 
recommend here the review article by Krieger [18] on the Wave Map problem which covers the most 
important results and references up to the year 2008. 

The main interest in studying the Wave Map equations and particularly the 2+1 Wave Map into a 2-sphere 
is the formation of singularities. These question is closely related to the previously discussed problem 
of global well-posedness or regularity. In the 1+1 case it is easy to show that no singularities can form, 
because the Wave Map equations are equivalent to the standard 1+1 wave equation. For spatial dimensions 
larger than two it is known that the Wave Maps can form singularities. See [ 1 8, 43] for detailed information 
on this. Only for the 2+1 case it was for a long time not known. In 2001 Bizori et. al. presented the 
first crucial numerical evidence for the singularity formation of the 2+1 dimensional Wave Map into the 
2-sphere. We are going to present their results and point out the important contributions during the last 
years in this exciting field. Based on their numerical observations of the equivariant case Bizoh et. al. 
formulated three conjectures about the singularity formation: 

Conjecture 1 (On blow-up for large data). For initial data (7) with sufficiently large energy, the solutions 
of equation (6) blow up infinite time in the sense that the derivative <9 r #(f,0) diverges as t /* T for some 
T>0. 

Conjecture 2 (On blow-up profile). Suppose that the solution $(?,r) of the initial-value problem (6) and 
(7) blows up at some time T > 0. Then, there exists a positive function s(t) \ Ofor t /*~T such that 

lim0(f,j(f)r) = # s (r) 

This conjecture was proven by Stmwe in [37], where, he not only proved the above statement, but further 
he showed that the existence of a non-trivial (non-constant) harmonic map is a necessary condition for a 
singularity formation. Struwe showed that the singularity formation appears as an energy concentration at 
the origin and the bubbling-qff of a harmonic map. The latter is the decomposition of the solution near the 
blow up time T in the following form (see [19]) 

(9) $(t,r) = V s ( S {t)r)+R(t,r). 

The function R(t,r) is a radiative term and #s (s(t)r) the rescaled static solution. The error term's local 
energy goes to zero as the time reaches the blow up time T. Equipped with the result in [ ] Krieger, 



2 

Small means here, small with respect to the energy. In other words, the energy is below a critical value. 
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Schlag and Tataru were able to construct the first blow up scenario for the 2+1 Wave Map into the 2-sphere 
[19]. This was later generalized to surfaces of revolution, which includes the 2-sphere, as target manifold. 

Conjecture 3 (On energy concentration). Suppose that the solution #(f,r) of the initial-value problem (6) 
and (7) blows up at some time T > 0. Define the kinetic and the potential energies at time t <T inside the 
past light-cone of the singularity by 



Ekin(t) 



T-t 



■& 2 rdr and E pot (t) = K 



T-t 



rdr 



Then: 



(a) the kinetic energy tends to zero at the singularity 

\imE kin (t)=0 

t/T 

(b) the potential energy equal to the energy of the static solution jJj concentrates at the singularity 

KmE pot (t)=E[# s )=4n. 



The behaviour of the energy, as described in the last conjecture is a direct consequence of the splitting (9) 
if one considers the vanishing energy of the radiation term in the past light cone. 

In a recent detailed study on the singularity formation Raphael and Rodnianski constructed a stable blow 
up scenario [26]. A different approach by Ovchinnikov and Sigal can be found in [ ]. In both publications 
the authors give an analytical form for the scaling function s(t). We are going to use this result later on 
in section 6.4. We want to mention here that the results by Raphael and Rodnianski cover all homotopy 
classes of the equivariant Wave Map. The singularity formation for higher homotopy classes was already 
studied by Rodnianski and Sterbenz in [30]. 

Extrinsic Formulation. In the later numerical simulations, we will use the extrinsic formulation, where 
we use Cartesian coordinates in the Euclidean space to avoid coordinate singularities. We write for the 
Wave Map: 

U : M 2+1 ->■ § 2 ^ R 3 

(U l (t,x,y),U 2 (t,x,y),U 3 (t,x,y)) . 

The line element simply is 

d4 3 = (df/ 1 ) 2 + (df/ 2 ) 2 + (df/ 3 ) 2 . 
The field equations, consisting of the equations of motion and the constraint equation, take the form: 

(10) U g U A ~2XU A =0, 

(11) U A U A -l=0. 

The Lagrangian parameter X =X[t,x,y) can be computed by multiplying the equation of motion ( 1 0) with 
ua and the use of the constraint equation (11) and its derivative. This leads to 

l = - l -{V"U B V a U B ) . 
The resulting non-linear wave equation has the constraint condition already imprinted 

a g u A + (v"u B v a u B ) u A = o. 

Now, we also introduce Cartesian coordinates on the base manifold M 2+1 



ds^i+i = d? 2 -dx 2 -dy 2 . 
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(16) X = ~\ 



By the relabelling u :— U , v := U 2 and w := t/ 3 we can write the equations (10) and (1 1) in the explicit 
form 

(12) ii~d xx u~dyyU — 2Xu = 1 

(13) v— dxxV— dyyV — 2Xv = 0, 

(14) W — dxxW — dyyW — 2Xw = 0, 

(15) <j){u,v,w) = u 2 + v 2 + w 2 -1=0. 
The Lagrangian multiplier is given explicitly by 

m 2 + v 2 + \v 2 — ((3 x m) 2 — (<3 x v) 2 — (d x w) 2 — (d y u) 2 — (dy.v) 2 — (<3yVv) / 
The hidden constraint for the velocities is the time derivative of (15): 

(17) \jf(u,v,w,u,v,w) = 2uu + 2vv + 2ww = 0. 

We use the following relation between the intrinsic and extrinsic formulation: 

(18) u = sin $ cos (p 

(19) v = sin#sin<p 

(20) w = cos & . 

Now, we can use these relations to write the static solutions (8) in the extrinsic formulation 

2x 2y 1 — r 2 

(21) Us (x >y ) = -—— > v s (x,y) = -—j, w$(x,y) = w s (r) = ±— -j, 

1 + r L 1 + r 1 1 + r z 

with the radial coordinate r = \] x 2 +y 2 . The ± sign in the expression for w$(x,y) corresponds to the ± 
sign in (8). It is easy to see that the static solutions, which are maps from R 2 into S 2 , are exactly the 
stereographic projections from the south (+) and the north (— ) poles, respectively. 

In conjecture 3 the existence of a scaling function s(t) is mentioned with the property to rescale the static 
solution in a way to approximate the dynamical solution around the origin. We can determine the scale 
factor at every fixed time as follows: we choose the component function w(t,r) of the Wave Map and 
choose s(t) in such a way that the second radial derivatives at the origin of the static and the rescaled 
dynamical solution (21) agree 1 . We use the function ws (r) with the + sign because we prescribe initial data 
in such a way that the stereographic projection from the south pole is singled out. The scaling relationship 

w s (r) = w(t,s(t)r) 

is assumed to be valid in a neighbourhood of the origin. Taking the second radial derivative and evaluating 
at the origin gives 

w^0)=s(t) 2 d rr w(t,0), 
from which we can determine the scaling factor as 

(22) s(t) 



V\drrW(t,0)\ ' 

4. Energy Conservation 

From (3) we can get the energy-momentum tensor by variation of the action with respect to the inverse 
metric g ab of the base space ^#: 

(23) Tab - VaU A V b U B 8 AB - ^VcU A V C U B 8 AB gab - tygab 

3 In [5] the first derivatives are used to compute the scaling function. But in the extrinsic formulation the function ws of the static 
solution is an even function and thus has vanishing odd derivatives at the origin. 
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In general, the energy of a matter field on a spatial slice in ^# as measured by an observer field with 
4- velocity t " is defined by 

E{t) = i f t a n b T ab sJ\y\A m x. 
2 Jl, 

The slicing of the manifold implies a (negative definite) spatial metric Yob on the s h ce which has the 
future -pointing normal vector field n b . For our purposes it is enough to consider flat space-like hyper-planes 
T, t . Furthermore, we choose static observers so that t" = (1,0,0) = n a . This and the energy-momentum 
tensor (23) lead to the following form of the energy contained inside a spatial domain QC1 2 

E$(t) = ^JjU A U A + d x U A d x U A + dyU A dyU A -X${U)] dxdy. 

Here we write E$ in order to indicate that we do not necessarily impose the constraint. This has the conse- 
quence that the Energy-Momentum tensor is not divergence-free, as we find by computing its divergence 

(24) V a T ah = -{d h l)^. 

From (24) we can see clearly that the energy-momentum tensor is divergence free, provided 0=0, i.e., 
if and only if the constraint is satisfied. But what will happen if this is not the case? This is what almost 
invariably happens in a numerical simulation. We can use the above formula to find the energy loss due to 
the constraint violation. 

We need the projection of (24) along the (constant) time-like vector t" 

(25) t b V>T ah = - (t b d b k) . 

We can interpret this equation as a continuity equation with an additional source on the right hand side. 
Using the fact that we are working on Minkowski space-time in Cartesian coordinates, we write for (25) 
with a dot as the time derivative 

(26) too + d% k = -k<j) 

where small Latin indices from the middle of the alphabet denote the spatial coordinates of the Minkowski 
space time. 

The required components of the energy-momentum tensor are 

Tm = \ U A U A - \ d k U A d k U A -A0=:p-A0 

(27) T Qk = U A d k U A =: j k 

with the 'true' energy density p of the system when the constraint is satisfied and the momentum current 
density j k . Now we can express (26) as 

p - A - X<j) + d k j k = — A 

p + d k j k = H 

and integrate over the domain £2 at a fixed time 

/ pdxdy+ d k j k dxdy= / X<j)dxdy. 
Jo. Ja Jn 

The first term on the left hand side in this equation is the time derivative of the 'true' energy E and the 
second term is an integral over the divergence of j k . With the help of Gauss' theorem, we can write the 
divergence term as a surface integral over the boundary <9£2 

E+ f n k j k dS= [ X<j>dxdy 
Jdn Ja 

(28) ==> E= [ Ai/Adxdy- / n k f dS 

Ja Jda 
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with the normal vector field n k on the boundary and the boundary element dS. The change of the energy 
between two instants of time is obtained by integrating over a time interval [t\,t2] 

(29) AE =E(t 2 ) -E{h) = [ 2 I lydtdxdy- [' [ n k j k dt dS . 

Jt { Jo. Jti Jdn 

The boundary term in (29) can be controlled by choosing appropriate boundary conditions. Inserting the 
definition of j k we find for the integrand 

(30) n k j k< ~=U A n k d k U A . 

In our numerical code we use homogeneous Neumann boundary conditions on all the unknowns Ua 

n k d k U A = . 

With these boundary conditions the boundary term (30) vanishes and (29) reduces to 

AE = [ f Xy/dt dxdy. 
Jtj Jn 

This is the change of the energy, which is due to the violation of the constraints. With other boundary 
conditions, for example Sommerfeld conditions, the boundary term would remain and describe an energy 
flux through the boundary. 



5. Numerical Setup 



5.1. Spatial Discretisation. In our investigations we use the Method of Lines approach for numerically 
solving PDEs. This means here we discretized the spatial variables with Finite Differences and evolve the 
resulting grid functions with either the fourth order Runge-Kutta or the Rattle method. 



(-1,1) (0,1) (1,1) 



1 
1 


£2 2 


(0,0) 




ill 





(-1-1) 



y 

A 



(1,0) 



(1,-1) 



FIGURE 1 . Domains of integration. Full domain £2j = [—1,1] x [—1,1] and the reduced 
(quarter) domain D.2 = [0, 1] x [0, 1]. 



In section 6.1, 6.2 and 6.3 we are going to compare the results for the time integration methods directly. 
There, we perform our simulations on the domain of integration Q.\ = [—1,1] x [—1,1] (see figure 1) and 
discretize it via an equidistant grid. 

x — >Xj = -\ + {j-\)h x j = l,...,N x 
y^y k = -l + (k-l)h y k = l,...,N y . 

with the grid spacings h x =Xj+\ —xj and h y =y k +\ —yk- All functions f(x,y) will be evaluated on the grid 
nodes: 



f(x,y) — > f(xj,yk)=fj.k- 
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This leads to approximations of the derivatives by differences of the function values at the grid points. For 
fourth order accuracy we write the following standard five point formulas for the derivative with respect to 
x at a point (j,k): 

~\ i*\ _ fj-2,k - + 8/,+l.j; - fj+2,k 

~j r I _ ~fj-2,k + 16/;- U - 30fjje + I6fj+l.k-fj+2.k 

wkm- 12h i x 

Derivatives with respect to y can be done equivalently. In all our computations we use the same grid 
spacing: h x = h y in x and y directions. 

For the simulations on the blow up dynamics in section 6.4, we use the symmetry of the system to reduce the 
domain of integration to a quarter of Q.\ only. This has the advantage, that we effectively double the number 
of grid points with the same numerical costs. The new grid covers therefore the domain Q.2 — [0, 1] x [0, 1] 
(see also figure 1). 

5.2. Boundary Conditions. As mentioned above, we use homogeneous Neumann conditions for the nu- 
merical calculations where we compare the Runge- Kutta and the Rattle method (this is important for the 
above described energy correction). We implement this boundary condition by imposing symmetries for 
the grid functions across the grid boundary. This gives us the possibility to determine the grid functions 
on points, which are beyond the boundaries (ghost points) which are needed for the evaluation of the finite 
difference operators. On the left boundary (j — 1) we get 

fo,k = f2,k and /_!,* = fi. k 

and on the right boundary (J = N x ) 

fN x +l,k = /v,-U an d fN x +2,k — fN x -2,k 

and similarly we set the boundary conditions in the y direction. 

The symmetries of the functions u and v make it necessary to impose different boundary conditions for the 
computations on the reduced domain £22 . For the function u, we need along the boundary x = and for 
the function v along y — Dirichlet boundary conditions and not Neumann ones. On those boundaries, 
we have u(0,y) — respectively v(x,Q) — 0. On the remaining three boundaries we can keep the Neu- 
mann conditions. To compute the ghost points for the Dirichlet boundary conditions we use the following 
relations on the left boundary 

fo,k = ~f2,k and /_!,*= -/ 3i * 

The boundaries where we have to change the boundary conditions are for both functions left boundaries, 
therefore, we not need the equivalent relations on the right boundaries. 

5.3. Initial data. The initial data we use are a polynomial bump with zeros at r\ and ri: 

, . [a f4 ( y )( V ) N )" forre[ri,r 2 l 
(31) M r ) = { ^ ( "~ n) > 

[0 otherwise. 

The constant A represents the amplitude of i9o(r). By using the relations (18), (19) and (20) we get the 
initial data for the extrinsic formulation. Throughout this article we use in all computations the parameters 
r\ = 0.5, r2 = 1 and n = 4. The initial data (31) describe a ring shaped bump in the xy-plane for the function 
w(t,x,y) centred around the origin. 

In order that the Cauchy problem is well posed we also need initial data for the velocities. We choose 

X 

u(0,x,y) = d r u = cos &o{r) #o(>") - 
v(0,x,y) = d r v = cos-dQ{r)$o(r)- r 
w(0,x,y) = d r w = — sm-&o(r) #o(r) • 
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With this choice of initial data for the velocities, the above described ring shrinks towards the origin. This 
means that the energy concentrates around the origin. After the function w(t,x,y) has reached a minimum 
close to the origin the ring starts expanding again. If the energy concentration at the origin is large enough, 
one expects a singularity formation like it is assumed in the three conjectures 1, 2, 3 above. 

Raphael and Rodnianski presented in [26] a set of initial data for analytical, but also numerical investi- 
gations. However we have chosen the initial data (31) to compare our numerical results with the ones in 
[5]. 

6. Results 

The main aim in this article is to compare the evolution of a constrained system by using an integration 
method for a free evolution and a method, which explicitly takes the constraint equations into account. As 
our model we choose the above described 2+1 dimensional Wave Map equations in the extrinsic formula- 
tion, as an example for a wave-like PDE system, which is constrained by additional (algebraic) conditions. 

In the case of the free evolution we use the equations (12), (13), (14) and replace the Lagrangian multiplier 
A with the expression (16). For this calculation we use the standard fourth order Runge-Kutta method 
(RK4). These results will be compared with the ones gained from the time evolution taking the constraint 
equation into account. Here we are going to use equations (12), (13), (14), (15) and (17). The time evolution 
will be done with the Rattle method (RTL) described in section 2. For the Rattle method we choose for the 
maximally allowed constraint violation the value 10~ 12 . This means that for all times t: ||0|| ma x( r ) < 10 -12 
as well as ||v / ||max(0 < 10~ 12 , where |j • |j max is the maximum norm. If this accuracy cannot be reached 
within 100 iterations the program aborts. For the computations presented here this was never the case. All 
following computations are done in the time interval t G [0, 1.6] and with a Courant-Friedrichs-Lewy factor 
At/h = 0.2. In the simulations, where we compare the two integration methods, we have chosen the values 
for the amplitude A in a way that none of the two integration methods breaks down. 

6. 1. Constraint functions. One of the main points in this article is the preservation of the constraint (15) 
and its derivative (17). In figure 2 we compare the constraint violation for the two different integration 
methods. We see clearly that the Runge-Kutta method leads to a permanent increase of the constraint 
violation during the time evolution. The Rattle method always keeps the constraint violation below the 
previously set value 10 -12 . 

For the velocity constraint yf we can see in figure 3 qualitatively the same results as for the constraint <p . 
But it is worth mentioning that the Rattle methods preserves the velocity constraint even a few orders better 
than requested. 

The accelerated increase of the constraint violation for the Runge-Kutta method at t sa 0.8 occurs at the 
moment when the wave packet is reflected at the origin. During this reflection process the spatial derivatives 
near the origin are very high. This also leads to a strong deviation of the function from its theoretically 
fixed value at the origin. 

We mention here, that the convergence for both methods is of the expected order. If we rescale the results 
shown in figure 2 and 3 for different numerical resolutions with the expected, fourth order, convergence 
rate, the curves for the Runge-Kutta method cannot visually be distinguished. To check the convergence 
of the Rattle method, we cannot use the constraint preservation, because in this method we predefine the 
maximum of the constraint violation. However, by checking the convergence rate with the function w with 
the function w we also obtain the expected second order convergence. 

6.2. Behaviour at the Origin. At the origin r = 0, respectively (x = 0,y = 0) the solution #(f , r) of the 
equivariant Wave Map equation (6) has the noteworthy feature that it is constant during the whole time 
evolution. One can show that from the equations that the solution evolving from our class of initial data 
necessarily satisfies for all t 

$(t,0) =0 and respectively w(t ,0,0) = 1 . 
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Figure 2. Maximum norm 1 10 1 |max(0 for the constraint during the time evolution for 
A = 0.4. 
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Figure 3 . Maximum norm 1 1 | m ax (0 for the velocity constraint \j/ during the time evo- 
lution for A = 0.4. The dotted line is the preset value for the maximum of the constraint 
violation. 

This condition is not enforced in our code. In figure 4 we compare how the two integration methods 
preserve this property. We can see that the Rattle method always keeps the deviation \w(t ,0,0) — 1 1 under 
the predefined maximum of 10~ 12 for the constraint violation. As expected from the previous results, 
the Runge-Kutta method cannot preserve the property, because it lets the solution drift away from the 
constraint manifold. In the literature, this drift is a well known phenomenon for constrained ODEs (see 
[21]). Numerical solutions of constrained systems, in the formulation where the Lagrangian multipliers are 
eliminated, show this behaviour in general. 

For larger values of A the Runge-Kutta method breaks down at the origin. This breakdown depends on 
the spatial resolution and is therefore a numerical phenomenon and not an indication for the collapse of 
the solution. The solution gained with the Rattle method shows a different behaviour. For large values 
of A the solution flips from w(f,0,0) = 1 to w(t ,0,0) = — 1. This indicates that it becomes increasingly 
difficult for the iterative algorithm of the constraint equation to find a solution and ultimately the iteration 
seems to approach the other static solution, the stereographic projection from the north pole, for which 
w(t, 0,0) = —1. We take this behaviour as an indication that the solution is in the blow-up regime. For 
the critical value of the amplitude, we take the value where the behaviour at the origin changes and find 
A* fa 0.81871173. 

6.3. Energy Conservation. An often used measure for the quality of numerical solutions of ODEs and 
PDEs is the energy conservation during the time evolution. Of course, this makes sense for non-dissipative 
systems only. We compare the energy conservation for the evolution computed with the Runge-Kutta and 
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FIGURE 4. Displacement |w(f,0,0) — 1| of the solution at the origin for A = 0.4 during 
the time evolution. The dotted line is the preset value for the maximum of the constraint 
violation. 



the ones obtained with the Rattle method. We can see in figure 5 that the average energy conservation with 
the Rattle method is better than with the Runge-Kutta method. 

In section 4 we presented the analytical formula for an energy loss, which is caused by the violation of the 
constraint equations. Our numerical investigations confirm the existence of this energy loss. We can use 
this term to correct the total energy E con = E — AE. From figure 5 we see that the corrected energy is better 
preserved than without correction, and also in average better than with the Rattle method. Only during the 
reflection at the origin the Rattle method results in an better energy conservation. 

The energy correction computed for the Rattle method is so small compared to the relative energy, that it 
can be neglected. 
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FIGURE 5. Relative energy |£ re i|(f) = |1 -E(t)/E(0)\ computed with the Runge-Kutta 
(corrected and uncorrected) and the Rattle methods for A = 0.4 and N = 641. 



6.4. Scaling Function. The conjectures by Bizon et.al. refer to the scaling function s(t), which we deter- 
mined in (22). This function plays a central role in the study of the singularity formation of Wave Maps 
[25, 26]. The importance comes from the fact, that with the help of s(t) the dynamical solution can be 
rescaled in a way that it approximates the static solution during the singularity formation, around the point 
where the singularity will appear. This means the scaling function is directly involved in the blow-up. The 
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fact that the scaling function plays an crucial role in the singularity formation was shown by Struwe in 
[37], The scaling function is an direct indicator for the shrinking of the static solution, which drives the 
singularity formation. 

In this part we determine the scaling function for different values of the criticality parameter A in its 
dependence on time. All the calculations in this subsection have been performed with the Rattle method 
because of its superior behaviour in the critical regime. As mentioned before the simulations for this part 
were computed on the reduced domain of integration £22. 

As we can see in figure 6, increasing the amplitude A towards the critical value A*, the time interval in 
which the scaling function is approximately constant, increases. This means, that during this period the 
solution remains almost constant, hovering in a quasi-static state. But for higher resolutions this states 
last shorter. Therefore we are very confident the method we presented is able to deal with such critical 
situations, but the need for much higher resolutions is evident. We interpret this quasi-static state as the 
fact that the numerical scheme cannot follow the shrinking of the harmonic map anymore. Therefore for 
this numerical resolution the critical situation occurred. When the scaling function increases again, the 
resolution is again sufficient to follow this process. 




Figure 6. The scaling function s(t) for different parameter values A. The value A* = 
0.81871173 is the critical amplitude. 



In figure 7 we can see the scaling function s(t) for different values of the amplitude parameter A. The black 
curve is the fit to the curve for the amplitude A = 0.81871 172, which is the one closest to the critical value 
A*. The interval t 6 [0.836,0.85] between the dotted lines is the domain used for the fitting procedure. The 
fit was done with respect to the analytic expression 

s(t) = ^-(T-t)exp(-y/-ln(T-t)+b) 

which is taken from [25]. The parameter b depends on the choice of the initial data. Our choice for the 
interval for the fit is a compromise between being close enough to the blow up time and the domain where 
we can be sure, that the scaling function numerically converges. The results presented here were computed 
with N = 641 grid points inx and y direction. For higher resolutions the amplitude A* = 0.81871 173 is not 
critical at all. The results of the fit were T = 0.94151363 and b = —2.02978334 with a residual error of 
1.31486579 • 10~ 8 . The fit parameter were computed with the Matlab function lsqcurvef it from the 
Optimization Toolbox. 

Figure 8 shows the rescaled function w(t,r/s(t)) at different times and the static solution w$(r). The 
amplitude is A = 0.81871172, very close to the critical value A* = 0.81871173. At t = 0.919375 the 
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Figure 7. Fit of the scaling function s(t) to the analytical expression, given in [25]. 
The computed blow up is T = 0.94151363. The value A = 0.81871 172 is the last under- 
critical value shown, i.e., for which the solution does not change its behaviour at the 
origin. 

solution reaches its overall minimum value of w = —0.91372693. Before this time the wave packet moves 
towards the origin, where it is reflected and subsequently moves away from the origin. We can see that the 
rescaled dynamic solution is approximated very well by the static solution near the origin. As mentioned 
before for larger amplitudes the solutions show a flip of the solution at the origin. In the caption of figure 
5 in [ ], Bizoh et. al. note the fact that the solution overshoots the value — %, which they consider as a 
necessary and sufficient condition for the blow up. We believe that in our case it is the observed flip at the 
origin, which indicates the same phenomenon. 



In this article we demonstrated that the Rattle method can be very useful for the numerical time evolution 
of relativistic field equations. In terms of constraint conservation it is, by construction, superior compared 
to a free evolution scheme like the standard fourth order Runge-Kutta method. The usage of the constraint 
equations also leads to a more accurate behaviour in extreme situations, like we showed it for the behaviour 
of the Wave Map solution at the origin. At this point the solution should always take the same value. The 
Rattle method preserves this symmetry in the accuracy previously set as the limit of the constraint violation. 
The Runge-Kutta method leads to a drift away from this state. 

We presented an energy correction term, which is a direct consequence of the constraint violation. The 
explicit computation of this deviation was used to correct the total energy. This energy deviation depends 
on the magnitude of the constraint violation, which is the reason why it can be ignored with the Rattle 
method by choosing the allowed constraint violation small enough. For the Runge-Kutta method it was 
possible to correct the energy in a way that the energy conservation during the whole time evolution is 
better than for the Rattle method. But this better energy conservation does not cure the other shortcomings 
namely the increasing constraint violation and violation of the symmetry at the origin. 

Finally we presented some results obtained in the equivariant case using the Rattle method. We computed 
the scaling function s(t) for different values for the amplitude and confirmed the expected behaviour that the 
rescaled dynamic solution approximates the static solution as the parameter A approaches a critical value. 
Since we obtain this behaviour without imposing the equivariance explicitly, this shows that the critical 
behaviour observed in the 1+1 codes is in fact stable under small non-spherical perturbations, which are 
inadvertently introduced due to the numerical errors. To our knowledge it has not been shown numerically 
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Figure 8. Ingoing wave packet for A = 0.81871172 rescaled with the scaling function 
s(t) for various values of t . The rescaled functions w(t,s(t)r) approximate the static solu- 
tion ws(r). The function w(t,r) reaches its minimum at t = 0.919375. The computation 
was done with N = 641 grid points on the quarter grid. 



that critical behaviour will also appear in the non-equivariant case. In a future article, we are going to show 
that under explicit violation of the equivariance, we still observe the critical behaviour. 

It is obvious from our results, that for an insight into the delicate details of the blow up dynamics much 
higher numerical resolutions are required. These could be achieved by using grid refinement techniques 
for example. However, the Rattle method applied to the extrinsic formulation of the 2+1 Wave Map seems 
to be a promising way to study this system without running into problems with coordinate singularities. 
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